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1 Introduction 



Sector decomposition in its practical aspect is a constructive method used to evaluate 
Feynman integrals numerically. The goal of sector decomposition is to decompose 
the initial integration domain into appropriate subdomains (sectors) and introduce, 
in each sector, new variables in such a way that the integrand factorizes, i.e. becomes 
equal to a monomial in new variables times a non-singular function. 

Originally it was used as a tool for analyzing the convergence and proving theorems 
on renormalization and asymptotic expansions of Feynman integrals El IH El E]- 
After a pioneering work [7] sector decomposition has become an efficient tool for 
numerical evaluating Feynman integrals (see Ref. [8] for a recent review). At present, 
there are two public codes performing the sector decomposition [9] and |10j . 

The latter one was named FIESTA which stands for "Feynman Integral Evaluation 
by a Sector decomposiTion Approach". Last year FIESTA has been greatly improved 
in various aspects [11] . The code is capable of evaluating many classes of integrals 
that one would not be able to evaluate with the original FIESTA 1. Moreover the 
code can now be applied to solve the problem of obtaining asymptotic expansions 
of Feynman integrals in various limits of momenta and masses and to find a list of 
all poles of an integral in space-time dimension d. During the last year FIESTA was 
widely used, some of application are listed in [T2] . 

In the current paper we present numerical results for master integrals (Mi's) for 
four-loop massless propagators which are relevant for many important physical ap- 
plications, like the calculation of the total cross-section of e + e~ annihilation into 
hadrons, the Higgs decay rate into hadrons, the semihadronic decay rate of the r lep- 
ton and the running of the fine structure coupling constant (see [131 E] for details). 
We confirm numerically the recent analytic results of work [TJ and evaluate an extra 
order in epsilon expansion for each ML 

2 Theoretical background and software structure 

FIESTA calculates Feynman integrals with the sector decomposition approach. It is 
based on the a-representation of Feynman integrals. After performing Dirac and 
Lorentz algebra one is left with a scalar dimensionally regularized Feynman integral 



where d = 4 — 2e is the space-time dimension, a n are indices, I is the number of 
loops and 1/E n are propagators. We work in Minkowski space where the standard 
propagators are the form l/(m 2 — p 2 — iO). Other propagators are permitted, for 
example, l/(v -k±iO) may appear in diagrams contributing to static quark potentials 



[15j 




(1) 



1 



or m hqet| where v is the quark velocity (see, e.g |16j ) . Substituting 



-i ain/2 roo 

1 daa^e-^, (2) 



E? T(a) 

changing the integration order, performing the integration over loop momenta, re- 
placing ati with Xit] and integrating over rj one arrives at the following formula (see 

e.g. my. 

F(a x , . . . ,a n ) = 

T(A-ld/2) f . f a,-A 



i=i 

where A = H" =1 a„ an d ^ an d -F are constructively defined polynomials of Xj. The 
formula fl3]) has no sense if some of the indices are non-positive integers, so in case of 
those the integration is performed according to the rule 

OO rf\ a ~ 1) 

^1vt/(*) = / (n) (0) 
o 1 (a) 

where a is a non-positive integer. 

After performing the decomposition of the integration region into the so-called 
primary sectors [7] and making a variable replacement, one results in a linear combi- 
nation of integrals of the following form: 

/ n' \ TjA-(l+l)d/2 

d Xi ... dx n , \ J] ] pA _ ld/2 (4) 



If the functions ^4x=757^~ had no singularities in e, one would be able to perform 
the expansion in e and perform the numerical integration afterwards. However, in 
general one has to resolve the singularities first, which is not possible for general U 
and F. Thus, one starts a process the sector decomposition aiming to end with a sum 
of similar expressions, but with new functions U and F which have no singularities 

(all the singularities are now due to the part IT7=i x>< j 3 )• Obviously it is a good idea 
to make the sector decomposition process constructive and to end with a minimally 
possible number of sectors. The way sector decomposition is performed is called a 
sector decomposition strategy and is an essential part of the algorithm. 

After performing the sector decomposition one can resolve the singularities by 
evaluating the first terms of the Taylor series: in those terms one integration is taken 
analytically, and the remainder has no singularities. Afterwards the e-expansion can 
be performed and finally one can do the numerical integration and return the result. 



3 Heavy-Quark Effective Theory 
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Please keep in mind that this approach works only using numerical integration: 
numeric values for all invariants should be specified at the very early stage, after 
generating the functions U and F. 

FIESTA is written in Mathematica [19] and C. The user is not supposed to use the 
C part directly as it is launched from Mathematica via the Mathlink protocol. When 
the integrand expressions are ready, Mathematica submits long strings representing 
integrands for integration; the C part translates them into an internal representation 
optimizing evaluation speed. Afterwards it uses some numerical integrator to perform 
the numerical integration of the integrand. The original FIESTA employed a Fortran 
implementation of Vegas as an integrator. Later we plugged in the Cuba library 
|18j . By default FIESTA uses the Vegas integrator, but this behavior can be easily 
controlled by the user. Both Mathematica and C parts can be efficiently parallelized 
on modern multi-core computers; the C part also parallelizable on clusters. 

The FIESTA user interface is based on Mathematica. To run FIESTA, the user has 
to load the FIESTA_2 . . . m into Mathematica 6 or 7. In order to evaluate a Feynman 
integral one has to use the command 

SDEvaluate [U , F , £ , indices , order] , 

where U and F are the functions from formula (j3J), £ is the number of loops, indices 
is the set of indices and order is the required order of the e-expansion. 

To avoid manual construction of U and F one can use a build-in function UF and 
launch the evaluation as follows: 

SDEvaluate [UF [loop_momenta, propagators , subst] , indices, order] , 

where subst is a set of substitutions for external momenta, masses and other values 
(please note that the code performs numerical integrations, therefore the functions U 
and F should not depend on any external kinematic invariants). 
Example: 

SDEvaluate [UF [{k} , {-k 2 , - (k+pi ) 2 , - (k+px+p 2 ) 2 , - (k+p!+p 2 +p 4 ) 2 } , 
{pi -^0,p 2 -^0,pl -)-0, pi p 2 -^-s/2,p 2 p 4 ->-t/2,pi p 4 -)-(s+t)/2, 
s-»-3,t-»-l}] , {1,1,1,1},0] 

performs an evaluation of the massless on-shell box diagram where the Mandelstam 
variables are equal to s = —3 and t = —1. 

3 Numerical results for four-loop massless propa- 
gators 

In [20] a full set of the four-loop massless propagator-like Mi's was identified. There 
exist 28 independent Mi's. Analytical results for these integrals were obtained in pQ. 
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M 36 , e l M 52 , e 1 M 43 , e 1 N ,e 2 

Figure 1: M 61 — M 43 : the thirteen complicated four- loop master integrals according 
to pQ. The integrals are ordered (if read from left to right and then from top to 
bottom) according to their complexity. The two Mi's M 52 and M 43 can be identically 
expressed through the three- loop nonplanar MI N . 



By an analytical result is meant not an analytical expression for a master integral 
taken at a generic value of the space-time dimension d (which is usually not possible 
except for the simplest cases), but rather the analytic expressions for proper number 
of terms in its Laurent expansion in d around the physical value d = 4. 

As it was shown in pQ , the full set of all 28 integrals can be divided in three parts: 

1. 9 "primitive" integrals which are expressible in terms of T-functions; 

2. 6 "simple" integrals which could be expressed in terms of T-f unctions and the 
so-called generalized two-loop diagram with insertions; 

3. remaining 13 "complicated" integrals are quite difficult for both analytical and 
numerical evaluation. 

The complicated Mi's are pictured in Fig. [TJ The Mi's are labeled as M y - where 
the first index, i, stands for the number of internal lines minus five while the second 
index, j, numerates (starting from 1) different integrals with the same i. e m after 
Mij stands for the maximal term in e-expansion of My which one needs to know for 
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evaluation of the contribution of the integral to the final result for a four-loop integral 
after reduction is done, see pp. That is, m stands for the maximal power of a spurious 
pole l/e m which could appear in front of My- in the process of reduction to masters. 

Primitive and simple integrals are known analytically. Two of the complicated 
integrals (M 43 and M 52 ) are related by a simple factor with the three-loop MI No P 
so it is enough to evaluate remaining eleven complicated Mi's M 61 — M 36 as well as 
first three terms of the ^-expansion of N Q . 

We have calculated all them by means of FIESTA. We used the Cuba Vegas in- 
tegrator with different parameters used for the numerical integration. Evaluations 
were performed on 8-core (2x4) Intel Xeon E5472 3.0 GHz, 4GB/core RAM, 4.6TB 
disk/node computers in fully parallel mode, i.e., both Mathematica and C parts were 
completely parallelized. The square of the external momentum q was chosen as -1: 
q 2 = —1. The FIESTA input for, say, the integral M 44 reads: 



Fm44= { 

-kl~2, -k3' 
-(kl-k2)-2, 

>; 



2, -k4~2, -(kl+q)~2, 
-(k2-k3)~2, -(k3-k4) 



(k2+q)~2, -(k4+q)~2, 



SDEvaluate[UF[{kl,k2,k3,k4}, Fm44,{q~2-> -1}], 
{1,1,1,1,1,1,1,1,1},!] ; 

Our results alongside with the corresponding analytical expressions (transformed 
to the numerical form) from [lj are presented in Tables 1 and 2. 

Within FIESTA it is implied that Feynman integrals are with the —k 2 — iO de- 
pendence of propagators and results are presented, in a Laurent expansion in e, by 
pulling out the factor m d / 2 e~ lE£ per loop, where 7^ is the Euler constant. Please, 
note that the overall normalization used by FIESTA is different from the one employed 
by the authors of pp. We denote by Mj a FIESTA result for an £-loop MI M,. The 
connection between both values reads: 



M 



Mi. 



Numerically, for 



Mi 



T(2 - 2e) 
4 the conversion factor is: 
8e + 36.710132 e 2 + 122.46185717e 3 + 329.99310668 e 4 + 



+ 758.778374e 5 + 1543.7276075 £ 6 + 2848.0962405 e 7 + 0(e 8 ) M, t (5) 



and for 



1 + 6 e + 21.53259889972766 e 2 + 0{e 3 



iV n . 



(6) 



5 



int. 


Degree 


Exact 


Cuba Vegas 


Cuba Vegas 


la. 


of s: 


Value: 


ouu uuu result. 


1 Kr\r\ nnn ^oonU- 
1 ouu uuu result. 




£ 4 


0.08333 


0.08333 ± 


0.08333 ± 






0.91666 


0.916668 ± 0.00003 


0.916667 ± 0.000018 




£ 2 


5.6425109 


5.64252 ± 0.00038 


5.64251 ± 0.00022 




£ 


Z ( .041z0ol 


Zrf.041o ± U.UUlo 


07 (3/1 1 4_ n nnn77 
ZI.041O it U.UUUrf I 


m 3 4 


£ 


QS fi^QOS 

yo.oo ( yzo 


qs r^s -i- n nnc;s 
yo.Doo it U.UUOo 


qs fi^s 4- n nn^/i 

yo.Doo it U.UUo4 




£ 




04Z. 1 00 It U.UZ1 


iAo I'xp. 4- n m 9 

O4Z.IO0 It U.U1Z 




£ 




fic;7 ss 4- n nsi 
oO 1 .00 it U.Uoi 


s K7 ss 4- n n/i s 
oO 1 .00 it U.U4o 




£ 3 


2659 825402 


2659 86 ± 32 


2659 84 ± 19 




5 4 




4344.27 ± 1.3 


4344.28 ± 0.75 




£ 


n fini nos 

U.OUlUZo 


n Rni n^n 4- n nnnno/i 

U.OU1UOU It U.UUUUZ4 


n fim nos 4- n nnnm 9 

U.OUlUZo It U.UUUU1Z 




£ 


( .4Z0U004 


7 /IOQO 4- n nnn/i 

1 .4ZOZ It U.UUU4 


7 /lOQi 4_ n nnno/i 

I.4Z01 It U.UUUZ4 


ivi 35 


£° 


44 91255 


44 9128 ± 0012 


44 9127 ± 00073 




£ X 


217.0209011 


217.026 ± 0.0062 


217.023 ± 0.0037 




£ 2 


780.4321125 


780.439 ± 0.022 


780.436 ± 0.013 




£ 3 




2678.18 ± 0.09 


2678.13 ± 0.053 




£- X 


5.1846388 


5.18467 ± 0.000072 


5.184645 ± 0.000042 


M 36 


£° 


38.8946741 


38.8950 ± 0.00068 


38.8948 ± 0.00039 




£ X 


240.0684359 


240.071 ± 0.0032 


240.069 ± 0.0019 




£ 2 




948.630 ± 0.016 


948.623 ± 0.0091 




s- 1 


20.7385551 


20.7386 ± 0.0004 


20.73860 ± 0.00023 


M41 


£° 


102.0326759 


102.034 ± 0.0051 


102.033 ± 0.003 




£ X 


761.5969858 


761.61 ± 0.019 


761.60 ± 0.011 




£ 2 




2326.21 ± 0.11 


2326.18 ± 0.062 



Table 1: Numerical results for the Mi's. In the third column the numerical values of 
the known analytical results are shown. The last two columns contain the results of 
evaluation on these integrals by FIESTA using the Cuba Vegas integrator with 500000 
and 1500000 sampling points correspondingly. For all Mi's we calculate one extra 
£-term (not known analytically). 
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Int. 


Degree 


Exact 


Cuba Vegas 


Cuba Vegas 


id: 


of e: 


Value: 


500 000 result: 


1 500 000 result: 




e- 1 


20.7385551 


20.7386 ± 0.00041 


20.73860 ± 0.00024 




e° 


145.3808999 


145.382 ± 0.0049 


145.381 ± 0.0029 




e 1 


985.9082306 


985.92 ± 0.023 


985.91 ± 0.014 




e 2 


— 


3930.68 ± 0.13 


3930.65 ± 0.076 


Maa 


e° 


55.5852539 


55.5858 ± 0.00054 


55.58537 ± 0.00031 




e 1 




175.325 ± 0.006 


175.325 ± 0.004 




e° 


52.0178687 


52.0184 ± 0.00052 


52.0181 ± 0.0003 


iK1 4o 


e 1 


175.496447 


175.50 ± 0.0062 


175.50 ± 0.0036 




e 2 




1475.29 ± 0.017 


1475.272 ± 0.0098 




P -\ 

c 


-o.ionjtjoo 


-O.lO^UIJ 31 U.UUUUol 


-o.ioiooi jz u.uuuu^o 


Mm 


e° 


-32.096143 


-32.0966 ± 0.00097 


-32.0962 ± 0.00057 




e l 


-91.1614758 


-91.157 ± 0.009 


-91.158 ± 0.0052 




e 2 




119.08 ± 0.074 


119.06 ± 0.043 




e° 


20.7385551 


20.7387 ± 0.00045 


20.73857 ± 0.00026 




e 1 


190.600238 


190.60 ± 0.004 


190.60 ± 0.0023 




e 2 


1049.194196 


1049.20 ± 0.025 


1049.20 ± 0.014 




e 3 




4423.86 ± 0.12 


4423.84 ± 0.072 




s' 1 


-10.3692776 


-10.36941 ± 0.00011 


-10.36931 ± 0.00006 




e° 


-70.99081719 


-70.989 ± 0.002 


-70.990 ± 0.0011 




£ l 


-21.663005 


-21.633 ± 0.023 


-21.650 ± 0.013 




e 2 




2832.86 ± 0.17 


2832.69 ± 0.096C1) 




e- 1 


-10.3692776 


-10.36940 ± 0.00011 


-10.36933 ± 0.00006 


M 62 


e° 


-58.6210462 


-58.6174 ± 0.0022 


-58.6187 ± 0.0013 




e 1 




244.69 ± 0.025 


244.681 ± 0.015 




e- 1 


-5.1846388 


-5.18470 ± 0.000078 


-5.18467 ± 0.000042 




e° 


14.397395 


14.40 ± 0.0014 


14.3989 ± 0.00081 




e 1 




740.00 ± 0.017 


739.979 ± 0.0099 



a Calculated with the Fortran Vegas using 1 550 000 samples. 

Table 2: Continuation of the table [U 
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Int. 


Degree 


Cuba Vegas 


Cuba Vegas 


id: 


of e: 


500 000 time: 


1 500 000 time: 




£- 4 


60.77s 


59.37s 




£" 3 


63.56s 


65.95s 




£" 2 


82.89s 


127.82s 




5- 1 


211.84s 


521.50s 


M 34 


£° 


401.53s 


1064.87s 




e 1 


586.88s 


1608.08s 




e 2 


988.27s 


2739.44s 




e 3 


1870.97s 


5225.72s 




e 4 


3422.80s 


9572.06s 




Total time: 


9847.39s 


23163.80s 




e~ 2 


64.29s 


76.14s 




e- 1 


178.16s 


426.75s 




e° 


325.72s 


855.19s 




e 1 


436.57s 


1169.30s 




e 2 


678.91s 


1828.51s 




e 3 


1275.64s 


3532.42s 




Total time: 


3764.31s 


8694.65s 




e" 1 


152.96s 


375.47s 




e° 


269.92s 


704.23s 




e 1 


354.42s 


959.94s 




e 2 


526.79s 


1442.97s 




Total time: 


1590.61s 


3769.25s 




s- 1 


185.32s 


274.82s 


M41 


e° 


691.48s 


1764.69s 




e 1 


928.23s 


2431.03s 




e 2 


1260.22s 


3379.72s 




Total time: 


3776.42s 


8562.06s 



Table 3: Timing for calculations of the Mi's. The last two columns contain time 
(in seconds) of numerical integration by the Cuba Vegas integrator with 500000 and 
1500000 sampling points. Also a total time for evaluation of each integral is given, 
including the Mathematica part. 
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Int. 


Degree 


Cuba Vegas 


Cuba Vegas 


id: 


01 £. 


oUU UUU time. 


1 ouu uuu time. 






176. 02 


246.99s 


M 42 


£ 


OoD.O / 


1 ( OZ.oUS 




£ 


91 f .90 


O/l Q K l 7K^ 
Z4oD. 1 OS 






1 oqq on 

IZOO.ZU 


Qosn o^ci 
ozo9.zos 




1 nt Q 1 tllTIf}" 

J-Uldl 11111c. 


^7^ Q9^ 


O'-tOO. VJ'-to 


7\ /f 

M44 


,-0 

£ 


7?no OHr, 

798. zUs 


2097.39s 




1 


1016.19s 


2713.11s 




lotal time: 


2174.60s 


5185.72s 




£° 


750.16s 


1906.93s 


M45 


e 1 


975.67s 


2533.32s 




2 

£^ 


1246.26s 


3256.41s 




Total time: 


3713.05s 


8416.63s 




1 

£ 1 


516.89s 


698.85 




n 

£° 


1676.80s 


4206.21 




1 

£ 


2881.73s 


7672.50 




£ 


O09 / .lOS 


9olo.los 




Total time: 


1U i 0O.U0S 


9/1977 Qnci 

Z4z / 1 .oUS 




£° 


42.13s 


104.29s 


iV 




75.01s 


201.78s 




2 

£ 


f\f~\ V T 

90.57s 


246.34s 




£ 


129. o4s 


Q /! 1 OOn 

o41.99s 




Total time: 


/ITT 1 On 

411.12s 


n^7 on^ 
9b 1 .2US 




£ -1 


2262.86s 


3495.28s 


M 61 


n 

£ u 


15242.10s 


39673.48s 




£ 


0l4ol.O0S 


1oz4oo.dzs 




£ 


2U2Ulo.olS 


1 <94o4U.UUs( j _J 




Total time: 


/ Oo/Z ^ .uus 






-l 

£ 


3003.05s 


4131.07s 


71 T 


^0 

£ 


16073.09s 


on/^nn pp^ 

39690.66s 




e 1 


63720.52s 


163026.12s 




Total time: 


156510.00s 


280778.00s 




e- 1 


273900.44s 


3316.17s 




£° 


14870.92s 


36434.93s 




e 1 


59206.88s 


151788.20s 




Total time: 


147870.00s 


262670.00s 



a Integration by the Fortran Vegas using 1 550 000 samples. 

Table 4: Continuation of the table [3j 
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A comparison with the analytical results shows that the integration with 500 000 
sampling points leads to the numerical result with 3-4 reliable digits in a quite rea- 
sonable time (see the tables [3] and HJ) while the integration with 1 500 000 sampling 
points reproduces the analytical results with 4-5 digits. We have also evaluated one 
extra term in the e-expansion of each MI which is currently unavailable analytically 
but is necessary for future five-loop calculations. For some technical reason^, for the 
highest e term of the integral M 6 i, we restricted ourselves with the value produced 
by the Fortran Vegas integrator which is not supported anymore. 

Surprisingly this planar integral (M 61 ) appears to be the most complicated one 
for numerical integration, see the table HI Non-planar integrals M 62 and M 63 are 
also complicated for FIESTA but much less than Mq\. Other integrals (including 
non-planars) are incomparably easier for numerical evaluation by FIESTA. 

The first thirteen Mi's from the Fig. [I] are very difficult for analytical evaluation, 
and only three of them had been checked in an independent way, see pQ. If even 
one of the remaining ten Mi's was evaluated incorrectly, it would change all physical 
results obtained with the use of these Mi's. 

Analytical results were obtained in [1] using so-called "glue-and-cut" symmetry to- 
gether with the procedure of reduction to Mi's. The reduction procedure is extremely 
complicated, it requires careful computer algebra programming and very large-scale 
computer evaluations. In the present paper we have performed the independent check 
of these Mi's using completely different approach, namely, sector decomposition, pro- 
viding a quite strong evidence for the correctness of the algorithms and their imple- 
mentation in pp. 

4 Conclusion 

Usually, analytical evaluation of multiloop MI is a kind of art. It requires a lot of 
efforts (and sometimes CPU time). In many situations, independent checkup is hardly 
any possible in reasonable time. That is why the simple in use tools for numerical 
evaluation like FIESTA are important. 

Some of the integrals presented in this paper are really complicated, and the orig- 
inal FIESTA 1 was not able to evaluate them at all. This was one of our motivations, 
in particular, to improve FIESTA so that it would cope with these (and, hopefully, 
many others) integrals. There had been both technical and theoretical complications 
which had to be solved [UJ for this aim. 

The successful check of the results of PP demonstrates that the current version of 
FIESTA is a powerful tool for evaluating integrals numerically and for cross-checking 
analytical results. 

The evaluation was performed before we've implemented the Cuba library. The integrator 
spends 1794640 seconds which is more than 20 days so we wouldn't like to load 8-core machine for 
such a period by the job which was already done 
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